program Integration;
{$APPTYPE CONSOLE}

uses
  SysUtils,
  Classes,
  Matrix in 'Matrix.pas',
  ODEModel in 'ODEModel.pas',
  ODEIntegrator in 'ODEIntegrator.pas',
  RungeKutta in 'RungeKutta.pas',
  TestModels in 'TestModels.pas',
  AdamsBashforthMoulton in 'AdamsBashforthMoulton.pas',
  Arenstorf in 'Arenstorf.pas';

{$R *.res}

var
  integrator : TODEIntegrator;
  model : TODEModel;

  input : char;

  InitialConditions : MatrixPtr;
  DensePoints : array of extended;
  TimeStart, TimeEnd, Step : extended;
  numberOfDensePoints : integer;
  i : integer;

  report1, report2, report3, report4 : TextFile;
begin
  { initializtion }

  writeln('1 for RK45 Arestorf; 2 for ABM Arenstorf; 3 for RK45 Bruccelator; 4 for ABM Bruccelator');
  readln(input);

  if( (input = '1') or (input = '2')) then
  begin
    initialConditions := CreateMatrix(4, 1);
    SetMatrixElement(initialConditions, 1, 1, 0.0);
    SetMatrixElement(initialConditions, 2, 1, -2.00158510637908252240537862224);
    SetMatrixElement(initialConditions, 3, 1, 0.994);
    SetMatrixElement(initialConditions, 4, 1, 0.0);

    TimeStart := 0.0;
    TimeEnd := 17.0652165601579625588917206249*2;
    Step := 1.0e-4;
  end;
  if( (input = '3') or (input = '4')) then
  begin
    initialConditions := CreateMatrix(2, 1);
    SetMatrixElement(initialConditions, 1, 1, 1.0);
    SetMatrixElement(initialConditions, 2, 1, 4.2665);
    TimeStart := 0.0;
    TimeEnd := 20.0;
    Step := 0.01;
  end;

  if(input = '1') then
  begin
    Step := 1.0e-2;

    numberOfDensePoints := Trunc((TimeEnd - TimeStart) / Step);
    SetLength(DensePoints, numberOfDensePoints);

    for i := 0 to numberOfDensePoints do
      DensePoints[i] := Step * i;

    model := TArenstorf.Create(initialConditions, 4);
    integrator := TRungeKutta45.Create(model, 1e-8, Step, TimeStart, TimeEnd, DensePoints);
  end;

  if(input = '2') then
  begin
    model := TArenstorf.Create(initialConditions, 4);
    integrator := TAdamsBashforthMoulton4.Create(model, 1e-3, 1e-3, TimeStart, TimeEnd);
  end;

  if(input = '3') then
  begin
    numberOfDensePoints := Trunc((TimeEnd - TimeStart) / Step);
    SetLength(DensePoints, numberOfDensePoints);

    for i := 0 to numberOfDensePoints do
      DensePoints[i] := Step * i;

    model := TBruccelator.Create(initialConditions, 2);
    integrator := TRungeKutta45.Create(model, 1.0e-12, 0.5, TimeStart, TimeEnd, DensePoints);
  end;
  if(input = '4') then
  begin
    Step := 1e-3;
    numberOfDensePoints := Trunc((TimeEnd - TimeStart) / Step);
    SetLength(DensePoints, numberOfDensePoints);

    for i := 0 to numberOfDensePoints do
      DensePoints[i] := Step * i;

    model := TBruccelator.Create(initialConditions, 2);
    integrator := TAdamsBashforthMoulton4.Create(model, 1e-12, Step, TimeStart, TimeEnd);
  end;

  integrator.Run;

  if(input='1') then begin
    AssignFile(report3, 'arenstorf-rk.txt');
    Rewrite(report3);
    for i := 0 to model.Y.Count-1 do
    begin
      writeln(report3,
            FloatToStr(GetMatrixElement(model.T[i], 1, 1)) + ' '
            + FloatToStr(GetMatrixElement(model.Y[i], 3, 1))
            + ' ' + FloatToStr(GetMatrixElement(model.Y[i], 4, 1))
            );
    end;
    CloseFile(report3);
  end;

  if(input = '2') then begin
    AssignFile(report4, 'arenstorf-abm.txt');
    Rewrite(report4);
    for i := 0 to model.Y.Count-1 do
    begin
      writeln(report4,
            FloatToStr(GetMatrixElement(model.T[i], 1, 1)) + ' '
            + FloatToStr(GetMatrixElement(model.Y[i], 3, 1))
            + ' ' + FloatToStr(GetMatrixElement(model.Y[i], 4, 1))
            );
    end;
    CloseFile(report4);
  end;

  if (input = '3' ) then
  begin
    AssignFile(report1, 'bruccelator-rk.txt');
    Rewrite(report1);

    for i := 0 to model.Y.Count-1 do
    begin
      writeln(report1,
            FloatToStr(GetMatrixElement(model.T[i], 1, 1)) + ' '
            + FloatToStr(GetMatrixElement(model.Y[i], 1, 1))
            + ' ' + FloatToStr(GetMatrixElement(model.Y[i], 2, 1))
            );
    end;
    CloseFile(report1);
  end;


  if (input = '4') then
  begin
    AssignFile(report2, 'bruccelator-abm.txt');
    Rewrite(report2);
    for i := 0 to model.Y.Count-1 do
    begin
      writeln(report2,
            FloatToStr(GetMatrixElement(model.T[i], 1, 1)) + ' '
            + FloatToStr(GetMatrixElement(model.Y[i], 1, 1))
            + ' ' + FloatToStr(GetMatrixElement(model.Y[i], 2, 1))
            );
    end;
    CloseFile(report2);
  end;
end.
